Ejemplo 2: Modelación de casos de COVID

Motivación

Debido a la pandemia producida por COVID-19, era de interés analizar posibles focos de infección y así poder tomar medidas para anticipar y evitar la propagación del virus.

En este caso mostraremos una posible alternativa de análisis para estimar riesgo de aparición del virus en diferentes radios censales de la Ciudad de Córdoba, Argentina.

Datos

Se trabajará con una base de datos vectorial (.gpkg) que incluye polígonos de los radios censales de Córdoba y Gran Córdoba. Para cada radio censal se cuenta con información de

  • Cantidad de habitantes en miles de personas (Poblacion)
  • Nivel de fragmentación urbano (fragmentacion)
  • Valor de la Tierra promedio
  • Porcentaje de hogares con nececsidades básicas insatisfechas (HogaresNBI)
  • Porcentaje de hogares con hacinamiento (HogariesHacinamiento)
  • Porcentaje de hogares con Jefe/Jefa con al menos nivel de estudio universitario (HogariesJefe.Univ.)
  • Porcentaje de hogares con computadoras (HogaresConComput.)
  • Cantidad de habitantes por km² (HabSup)
  • Cantidad de bancos dentro del radio censal (Bancos)
  • Cantidad de farmacias dentro del radio censal (Farmacias)
  • Cantidad de centros de salud dentro del radio censal (EstSaludInt)
  • Cantidad de bancos, farmacias y centros de salud del radio censal (Ban_Farm_Salud)
  • Cantidad de casos detectados de COVID-19 (Casos)

Lectura base de datos

Carga paquetes
library(sf)
library(spdep)
library(tmap)
library(INLA)

# library(mapview)
# library(leaflet)
# library(leafem)
# library(sp)
#
# 
# library(ggplot2)
# library(ggpubr)
# library(knitr)
datos <- st_read("data/Base_07_07_radios.gpkg", quiet = TRUE)
  • Cálculo de número de casos esperados por radio censal
datos$E <- datos$Poblacion * sum(datos$Casos) / sum(datos$Poblacion)
tm_shape(datos) +
  tm_polygons(col = 'E')

  • Cálculo de cociente de infección estandarizada (Standardized Infection Ratio, SIR) por radio censal
datos$SIR <- datos$Casos / datos$E
tm_shape(datos) +
  tm_polygons(col = 'SIR')

  • Cálculo de valor de la tierra relativo al promedio
datos$ValorTierra <- datos$ValorTierra / mean(datos$ValorTierra)
  • Referencia para el efecto aleatorio y el término del error
datos$re_u <- 1:nrow(datos)
datos$re_v <- 1:nrow(datos)

Ajuste del Modelo

Generación de lista con vecindarios

# Generación de lista con vecindarios
nb <- poly2nb(datos)

head(nb)

plot(st_geometry(datos), border = "grey", lwd = 0.5)
plot(
  nb,
  coordinates(as(datos, "Spatial")),
  add = TRUE,
  col = "blue",
  points = FALSE,
  lwd = 0.5
)

## [[1]]
##  [1]    2    3    4   29   37  237  251  252  258  266  267 1498 1500
## 
## [[2]]
## [1]  1  3  8 14 37
## 
## [[3]]
## [1] 1 2 4 8
## 
## [[4]]
## [1]   1   3   5   6   7   8 266 267
## 
## [[5]]
## [1]   4   6 267
## 
## [[6]]
## [1]   4   5   7  10  11 270
# Generación de vecindarios para INLA
file_adj <- tempfile("map.adj1")
nb2INLA(file_adj, nb)
g <- inla.read.graph(filename = file_adj)
# Ajuste del Modelo inflado en ceros
summary(datos)


formula  <-
  Casos ~ 
  Fragmentacion +  
  HogaresNBI + 
  HogaresHacinamiento +
  HogaresJefe.Univ. + 
  Bancos + 
  f(re_u, model = "besag", graph = g) + 
  f(re_v, model = "iid")

# formula  <-
#   Casos ~ Bancos + ValorTierra +  f(re_u, model = "besag", graph = g) + f(re_v, model = "iid")

res_inflpoi <-
  inla(
    formula,
    family = "zeroinflatedpoisson1",
    data = datos,
    E = E,
    control.predictor = list(compute = TRUE)
  )
summary(res_inflpoi)
##    Poblacion      Fragmentacion    ValorTierra         HogaresNBI     
##  Min.   :   7.0   Min.   :1.000   Min.   : 0.01171   Min.   : 0.0000  
##  1st Qu.: 629.8   1st Qu.:2.000   1st Qu.: 0.14349   1st Qu.: 0.3755  
##  Median : 835.0   Median :4.000   Median : 0.35141   Median : 1.0855  
##  Mean   : 873.2   Mean   :3.367   Mean   : 1.00000   Mean   : 1.6939  
##  3rd Qu.:1086.0   3rd Qu.:4.000   3rd Qu.: 0.73210   3rd Qu.: 2.4040  
##  Max.   :3273.0   Max.   :4.000   Max.   :12.88503   Max.   :20.6897  
##  HogaresHacinamiento HogaresJefe.Univ. HogaresConComput.     HabSup       
##  Min.   :0.0000      Min.   : 0.0000   Min.   : 0.00     Min.   :  0.006  
##  1st Qu.:0.1364      1st Qu.: 0.9802   1st Qu.:11.96     1st Qu.: 44.619  
##  Median :0.5051      Median : 3.3943   Median :18.05     Median : 73.769  
##  Mean   :0.8617      Mean   : 5.5653   Mean   :20.02     Mean   : 83.815  
##  3rd Qu.:1.2639      3rd Qu.: 9.0126   3rd Qu.:25.26     3rd Qu.: 97.168  
##  Max.   :6.6606      Max.   :26.6055   Max.   :57.17     Max.   :606.024  
##      Bancos          Farmacias       EstSaludInt     Ban_Farm_Salud  
##  Min.   :0.00000   Min.   :0.0000   Min.   :0.0000   Min.   :0.0000  
##  1st Qu.:0.00000   1st Qu.:0.0000   1st Qu.:0.0000   1st Qu.:0.0000  
##  Median :0.00000   Median :0.0000   Median :0.0000   Median :0.0000  
##  Mean   :0.05301   Mean   :0.2446   Mean   :0.1536   Mean   :0.4512  
##  3rd Qu.:0.00000   3rd Qu.:0.0000   3rd Qu.:0.0000   3rd Qu.:1.0000  
##  Max.   :8.00000   Max.   :6.0000   Max.   :6.0000   Max.   :9.0000  
##      Casos                    geom            E                 SIR        
##  Min.   : 0.0000   MULTIPOLYGON :1660   Min.   :0.002241   Min.   : 0.000  
##  1st Qu.: 0.0000   epsg:22174   :   0   1st Qu.:0.201595   1st Qu.: 0.000  
##  Median : 0.0000   +proj=tmer...:   0   Median :0.267299   Median : 0.000  
##  Mean   : 0.2795                        Mean   :0.279518   Mean   : 1.007  
##  3rd Qu.: 0.0000                        3rd Qu.:0.347649   3rd Qu.: 0.000  
##  Max.   :36.0000                        Max.   :1.047749   Max.   :89.967  
##       re_u             re_v       
##  Min.   :   1.0   Min.   :   1.0  
##  1st Qu.: 415.8   1st Qu.: 415.8  
##  Median : 830.5   Median : 830.5  
##  Mean   : 830.5   Mean   : 830.5  
##  3rd Qu.:1245.2   3rd Qu.:1245.2  
##  Max.   :1660.0   Max.   :1660.0  
## 
## Call:
##    c("inla.core(formula = formula, family = family, contrasts = contrasts, 
##    ", " data = data, quantiles = quantiles, E = E, offset = offset, ", " 
##    scale = scale, weights = weights, Ntrials = Ntrials, strata = strata, 
##    ", " lp.scale = lp.scale, link.covariates = link.covariates, verbose = 
##    verbose, ", " lincomb = lincomb, selection = selection, control.compute 
##    = control.compute, ", " control.predictor = control.predictor, 
##    control.family = control.family, ", " control.inla = control.inla, 
##    control.fixed = control.fixed, ", " control.mode = control.mode, 
##    control.expert = control.expert, ", " control.hazard = control.hazard, 
##    control.lincomb = control.lincomb, ", " control.update = 
##    control.update, control.lp.scale = control.lp.scale, ", " 
##    control.pardiso = control.pardiso, only.hyperparam = only.hyperparam, 
##    ", " inla.call = inla.call, inla.arg = inla.arg, num.threads = 
##    num.threads, ", " blas.num.threads = blas.num.threads, keep = keep, 
##    working.directory = working.directory, ", " silent = silent, inla.mode 
##    = inla.mode, safe = FALSE, debug = debug, ", " .parent.frame = 
##    .parent.frame)") 
## Time used:
##     Pre = 1.49, Running = 21.4, Post = 0.657, Total = 23.6 
## Fixed effects:
##                       mean    sd 0.025quant 0.5quant 0.975quant   mode kld
## (Intercept)         -0.662 0.492     -1.641   -0.658      0.287 -0.648   0
## Fragmentacion        0.169 0.106     -0.036    0.168      0.378  0.167   0
## HogaresNBI          -0.052 0.070     -0.199   -0.050      0.078 -0.044   0
## HogaresHacinamiento  0.115 0.140     -0.154    0.113      0.397  0.108   0
## HogaresJefe.Univ.    0.052 0.023      0.008    0.052      0.097  0.052   0
## Bancos               0.138 0.194     -0.249    0.140      0.514  0.145   0
## 
## Random effects:
##   Name     Model
##     re_u Besags ICAR model
##    re_v IID model
## 
## Model hyperparameters:
##                                                          mean     sd 0.025quant
## zero-probability parameter for zero-inflated poisson_1  0.598  0.052      0.503
## Precision for re_u                                      0.346  0.062      0.242
## Precision for re_v                                     31.850 40.820      5.001
##                                                        0.5quant 0.975quant
## zero-probability parameter for zero-inflated poisson_1    0.595      0.706
## Precision for re_u                                        0.340      0.485
## Precision for re_v                                       19.933    134.836
##                                                          mode
## zero-probability parameter for zero-inflated poisson_1  0.581
## Precision for re_u                                      0.327
## Precision for re_v                                     10.215
## 
## Marginal log-Likelihood:  -2315.06 
##  is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')

Visualización Espacial

Visualización Interactiva

datos$RR <- res_inflpoi$summary.fitted.values[, "mean"]
datos$RR_LI <- res_inflpoi$summary.fitted.values[, "0.025quant"]
datos$RR_LS <- res_inflpoi$summary.fitted.values[, "0.975quant"]

popup_vars <- c(
  "Fragmentación (Nivel)" = 'Fragmentacion',
  "Valor Tierra Medio ($/m2)" = 'ValorTierra',
  "Hogares NBI (%)" = 'HogaresNBI',
  "Hogares Hacinamiento (%)" = 'HogaresHacinamiento',
  "Hogares Jefe Univ.(%)" = 'HogaresJefe.Univ.',
  "Bancos",
  "Casos",
  "Población" = 'Poblacion',
  "E",
  "SIR",
  "RR_LI",
  "RR",
  "RR_LS"
)
tmap_mode('view')
mapas <-
  tm_basemap(c(Urbano = "OpenStreetMap", Satelite = "Esri.WorldImagery")) +
  tm_shape(datos,
           name = 'Casos')  +
  tm_polygons(
    col = "Casos",
    id = "Casos",
    border.col = "gray50",
    border.alpha = .5,
    style = "fixed",
    title = "Casos por Radio",
    palette = "YlOrBr",
    breaks = c(0, 1, 2, 3, 5, 10, 15, 25, 36),
    popup.vars = popup_vars,
    legend.format = list(
      scientific = TRUE,
      format = "f",
      digits = 0
    )
  ) +
  
  tm_shape(datos,
           name = 'Tasa de Infección - SIR')  +
  tm_polygons(
    col = "SIR",
    id = "SIR",
    border.col = "gray50",
    border.alpha = .5,
    style = "fixed",
    palette = hcl.colors(7, "ag_GrnYl"),
    breaks = c(0, 5, 10, 15, 20, 25, 30, 50, 60, 90),
    popup.vars = popup_vars,
    legend.format = list(
      scientific = TRUE,
      format = "f",
      digits = 0
    )
  ) +
  
  tm_shape(datos,
           name = 'Riesgo Relativo')  +
  tm_polygons(
    col = "RR",
    id = "RR",
    style = "fixed",
    palette = "viridis",
    alpha = 0.8,
    title = "Riesgo Relativo",
    breaks = c(0.4, 1, 2, 4, 8, 10, 15, 20, 30, 50, 85, 96),
    popup.vars = popup_vars,
    border.col = "gray50",
    border.alpha = .5,
    legend.format = list(
      scientific = TRUE,
      format = "f",
      digits = 1
    )
  )
mapas
img <-
  "https://ig.conae.unc.edu.ar/wp-content/uploads/sites/68/2022/04/cropped-cromas-68-1.png"

map <-
  tmap_leaflet(mapas) %>%
  leafem::addLogo(img,
                  url = "https://ig.conae.unc.edu.ar/",
                  width = 136,
                  height = 50) %>%
  leaflet::addMiniMap(position = "bottomleft",
                      width = 150,
                      height = 150)
map

# mapshot(map, paste0("Mapa_radios_", format(Sys.time(), "%d_%m_%Y"), ".html"))

Visualización Estática

tmap_mode('plot')

boundary_box <- matrix(c(4357000, 6508000, 4405900, 6560000))

mapacasos <-
  tm_shape(datos, bbox = boundary_box)  +
  tm_polygons(
    col = "Casos",
    id = "Casos",
    border.col = "gray50",
    border.alpha = .5,
    style = "fixed",
    title = "Casos por Radio",
    palette = "YlOrBr",
    breaks = c(0, 1, 2, 3, 5, 10, 15, 25, 36),
    legend.format = list(
      scientific = TRUE,
      format = "f",
      digits = 0
    )
  ) + 
  tm_scale_bar(position = c("right", "bottom")) +
  tm_compass(type = "8star", position = c("left", "bottom"))

mapacasos

mapaSIR <-
  tm_shape(datos, bbox = boundary_box)  +
  tm_polygons(
    col = "SIR",
    id = "SIR",
    border.col = "gray50",
    border.alpha = .5,
    style = "fixed",
    palette = "YlOrBr",
    breaks = c(0, 5, 10, 15, 20, 25, 30, 50, 60, 90),
    legend.format = list(
      scientific = TRUE,
      format = "f",
      digits = 0
    )
  ) + 
  tm_scale_bar(position = c("right", "bottom")) +
  tm_compass(type = "8star", position = c("left", "bottom"))

mapaSIR

mapaRR <-
  tm_shape(datos, bbox = boundary_box)  +
  tm_polygons(
    col = "RR",
    id = "RR",
    style = "fixed",
    palette = "YlOrBr",
    title = "Riesgo Relativo",
    breaks = c(0.4, 1, 2, 4, 8, 10, 15, 20, 30, 50, 85, 96),
    border.col = "gray50",
    border.alpha = .5,
    legend.format = list(
      scientific = TRUE,
      format = "f",
      digits = 1
    )
  ) + 
  tm_scale_bar(position = c("right", "bottom")) +
  tm_compass(type = "8star", position = c("left", "bottom"))

mapaRR